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ABSTRACT 


A  molecular  dynamics  computer  simulation  was  used  to  investigate  several 
techniques  of  generating  liquid  Cu  targets.  The  target  with  the  best  liquid 
characteristics  was  subjected  to  one  KeV.  argon  ion  bombardment  as  a  preliminary 
study  of  the  sputtering  of  liquids.  The  techniques  of  warming  by  impulse  and 
warming  by  initially  displacing  atoms  from  their  equilibrium  positions  were 
compared.  Both  methods  produced  targets  with  good  liquid  properties.  The  energy 
became  equally  partitioned  between  kinetic  and  potential  energy  and  all  targets 
equilibrated  within  400  fs.  The  range  of  a  typical  atom  during  the  time  of 
equilibration  was  found  to  be  restricted  to  its  initial  neighborhood.  The  preliminary 
sputtering  study  resulted  in  a  sputtering  yield  increase  of  40%  over  the  solid  target, 
for  a  low  index  crystal  plane. 
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I.  INTRODUCTION 


A.  HISTORICAL  OVERVIEW  OF  SPUTTERING 

In  1851,  Plucker  [ref.  1]  observed  that  the  gas  in  x— rays  tubes  was  continually 
removed.  He  attributed  this  phenomenon  to  the  ionization  of  residual  gases  and 
their  consequent  absorption  by  the  inner  surface  of  the  x— ray  tube.  In  1852,  Grove 
[ref.  2]  noticed  that  the  x— ray  tube  surface  struck  by  these  particle  ions  was  eroded 
due  to  the  removal  of  target  material.  He  called  this  phenomena  cathode  sputtering. 

At  that  point  in  history,  the  action  of  target  erosion  or  sputtering  was 
considered  undesirable  in  lab  equipment,  and  was  only  interesting  so  long  as  one 
could  minimize  its  unwanted  effects.  Consequently,  sputtering  was  not  investigated 
systematically  for  more  than  fifty  years. 

About  fifty  years  after  Grove  reported  his  cathode  sputtering  findings. 
Goldstein  [ref.  3]  presented  conclusive  evidence  that  sputtering  was  indeed  caused 
by  the  positive  atoms  of  the  discharge  impacting  on  the  cathode  target. 

In  190S.  Stark  [ref.  4]  advanced  the  concept  of  the  individual  sputtering  event 
on  an  atomic  scale.  He  developed  a  collision  theory  which  treated  sputtering  as  a 
sequence  of  binary  collisions  initiated  by  the  collision  of  one  incident  ion  at  a  time. 
He  also  presented  a  second  theory  called  the  hot-spot  model  which  attributed  the 
sputtering  action  to  localized  high  temperature  heating  of  the  target  and 
evaporation  of  atoms.  In  1921,  Thompson  [ref.  5]  proposed  that  sputtering  was 
caused  by  the  release  of  radiation  as  the  bombarding  ions  struck  the  target.  The 
following  year  Bush  and  Smith  [ref.6]  suggested  that  sputtering  was  caused  by  the 
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expansion  of  gas  adsorbed  by  the  target  material,  and  Kingdon  and  Langmuir 
[refs.  7.8]  conducted  a  sputtering  experiment  which  yielded  ejecta  from  the  surface 
layers  of  the  target.  They  bombarded  thoriated  tungsten  with  ions  in  a  glow 
discharge  tube.  The  experiment  demonstrated  that  the  bulk  of  the  ejecta  was 
coming  from  the  thin  film  of  thorium  on  the  target's  surface  instead  of  the  tungsten 
substrate.  Their  results  suggested  a  momentum  transfer  ejection  mechanism  for 
sputtering. 

In  1926.  Von  Hippie  and  Blechschmidt  [refs.  9— 12J.  proposed  a  theory  that 
described  sputtering  as  an  evaporation  of  the  surface  atoms.  Von  Hippie  extended 
Stark's  hot-spot  model  and  attempted  to  formulate  a  sputtering  theory  on  the  basis 
of  local  heating.  He  expressed  the  view  that  local  heating  was  the  only  feasible  way 
to  explicitly  treat  the  statistics  of  the  complex  collisions  occurring  in  a  sputtering 
event. 

Approximately  ten  years  later.  Lamar  and  Compton  [ref.  13]  published  .4 
Special  Theory  of  Cathode  Spattering  which  led  to  the  thermal  spike  concept.  The 
thermal  spike  was  based  on  momentum  transfer  between  the  incomming  ion  and  the 
lattice  atoms.  This  theory  suggested  that  a  long-lived  high  temperature  volume 
persisted  in  the  target  even  after  the  collision  cascade  was  completed. 

In  1931,  Guntherschulze  and  Meyer  [refs.  14,15]  were  the  first  to  recognized  the 
importance  of  minimizing  ambient  pressure  in  sputtering  experiments.  They  used  a 
high  vacuum  tube  and  took  the  precaution  of  removing  several  top  layers  from  the 
target  to  clean  its  surface.  Consequently,  their  results  were  the  first  to  satisfy  the 
conditions  for  reproducible  sputtering  yield  determination.  About  ten  years  later. 
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Penning  and  Mobious  [ref.  Hi]  conclusively  demonstrated  the  significant  effects  of 
ambient  pressure  on  the  sputtering  yield. 

By  the  early  1950's,  a  renewed  interest  in  sputtering  coalesced  in  the  scientific 
community.  This  was  brought  about  primarily  by  the  fact,  that  sputtering  was 
found  to  have  technological  value.  It  was  discovered  that  the  ion— trapping  process 
in  sputtering  could  be  used  as  a  pumping  effect  for  low  pressure  electronic  systems. 
It  was  also  found  that  sputtering  could  be  used  to  clean  target  surfaces.  Sputtering 
became  useful  to  industry  and  its  status  was  raised  from  that  of  a  laboratory 
hindrance  to  one  meriting  serious  scientific  research  effort. 

In  1952.  Keywcll  [ref.  17]  formulated  a  sputtering  theory  which  made  use  of 
existing  neutron  transport  theories  originally  developed  for  nuclear  reactor  design. 
Key  well's  work,  as  well  as  subsequent  calculations  by  Harrison  [ref.  IS]  made 
important  contributions  to  the  field  of  sputtering  by  introducing  probability 
concepts  in  terms  of  coll isiotial  cross  sections. 

In  1951.  Wehner  [ref.  19]  published  his  findings  on  crystal  structure  effects  in 
the  flux  of  sputtered  particles.  This  was  a  significant  advance  in  that  it  proved 
conclusively  that  Stark’s  hot-spot  model  was  incapable  of  fuily  explaining  the 
sputtering  mechanism.  Wehner's  conclusions  fortified  collisional  theory  as  an 
important  part  of  sputtering  theory.  His  findings  created  great  interest  in  ihe 
collisional  aspect  of  sputtering  and  served  to  point  the  way  for  following  research 
efforts. 

In  the  decade  that  followed,  (mid  50's  to  mid  60's)  the  main  emphasis  in 
sputtering  research  was  directed  towards  the  study  of  crystal  lattice  effects.  Such 
studies  produced  two  main  theories.  These  were  the  channeling  theory  and  the 


focusing  collision  theory.  The  channeling  theory  was  successful  in  describing  the 
angular  variation  of  the  yield  when  the  incident  ions  where  closely  aligned  to  one  of 
the  target's  principal  axis  or  planes.  The  theory  failed  for  general  alignment  of  the 
incident  ions.  The  focusing  collision  theory  was  the  product  of  Silsbee  [ref.  20].  His 
theory  showed  that  momentum  was  transferred  within  the  target's  crystal  structure 
along  preferred  directions.  He  introduced  the  idea  of  momentum  transport  without 
requiring  mass  displacement.  This  concept  was  later  refered  to  as  focusons. 
Focusing  collision  theory  was  relatively  successful  in  modelling  the  observed  ejection 
patterns  from  materials  with  a  high  degree  of  symmetry,  such  as  cubic  crystals,  but 
only  if  bombarded  with  high  energies.  This  theory,  however,  failed  to  match  the 
W’elmer  spots  and  proved  to  be  particularly  inadequate  for  targets  with  low 
symmetry  such  as  hexagonal  crystal  structures.  Experimental  results  demonstrated 
that  the  sputtering  yield  was  significantly  dependent  on  the  crystallographic 
orientation  of  the  target  with  respect  to  the  incident  ion  beam. 

The  theory  of  collisional  cascades  was  further  advanced  by  the  work  of  Sigmund 
(refs.  21,22].  He  used  Lindhard's  range  theory  [ref.  23]  and  proposed  that  for 
amorphous  solids,  the  collision  cascades  could  be  described  by  a  Boltzmann 
•transport  equation  [ref.  24],  His  equations  required  the  target  surface  to  have  a 
disordered  structure,  not  the  long  straight  rows  of  atoms  intersecting  the  surface 
used  in  the  focusing  model.  He  obtained  first  order  asymptotic  solutions  for  cascade 
events  produced  by  ions  of  medium  to  large  atomic  number  with  energy  in  the  KeV. 
His  theory  became  the  reference  standard  for  sputtering  yield  measurements  despite 
its  limitation  of  application  to  polycrystalline  solids  with  randomly  oriented 
crystallites  as  an  amorphous  approximation.  This  idea  was  simultaneously 


investigated  by  Thompson  [ref.  25).  He  proposed  that  the  ejected  atom  was  affected 
by  the  surface  attraction  of  the  target's  surface,  causing  a  deviation  on  the  velocity 
vector  of  the  outgoing  particle.  This  resulted  in  a  distortion  of  the  angular 
distribution  of  ejected  particles.  Another  important  idea  developed  from  collisional 
cascades  was  that  of  replacement  collision  sequences.  In  replacement  collision 
sequences,  the  moving  atom  replaces  an  atom  on  its  lattice  site.  The  vacated  atom 
the  proceeds  to  strike  and  replace  the  next  atom  in  the  row.  The  sequence 
propagates  as  each  atom  replaces  the  next. 

Despite  its  long  history,  sputtering  researchers  have  not  been  able  to  formulate 
an  analytical  theory  which  can  fully  explain  sputtering.  There  are  still  differences 
between  the  theoretical  predictions  and  experimental  data,  and  it  seems  that  an 
analytically  simple  form,  in  all  likelyhood,  will  not  be  sufficiently  flexible  to 
correctly  describe  all  aspects  of  sputtering.  The  crux  of  the  problem  lies  in  the  fac* 
that  in  many  cases,  the  theory  requires  the  solution  of  a  many-body  problem  with 
multiple  interactions.  This  is  a  formidable  task  and  one  which  historically  has 
yielded,  at  best,  only  approximate  solutions.  Another  approach  was  required  which 
could  explain  and  predict  sputtering  events  based  on  a  few  simple  laws.  With  the 
advent  of  the  high  speed  coputer  capable  of  handling  the  voluminous  amount  of 
required  calculations,  some  researchers  turned  to  computer  simulations. 

B.  COMPUTER  SIMULATIONS  OF  SPUTTERING 

I.  Historical  Overview 

By  the  late  i960's,  the  once  curious  laboratory  observation  of  Plucker 
[ref.  1]  and  Grove  [ref.  2],  had  been  studied  for  approximately  100  years  and 
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sputtering  theories  had  evolved  into  dependable  sputtering  yield  predictors. 
Especially  successful  in  determining  sputtering  yield  were  the  statistical  theories  of 
Thompson  [ref.  25],  ejected  atom  energy  distribution  function  dN/dE,  and  Sigmund 
[ref.  24],  sputtering  yield  as  a  function  of  energy,  Y(E). 

The  advent  of  the  high  speed  computer  made  it  practical  to  perform  the 
many  calculations  required  to  treat  the  collision  cascade  problem  in  a  three 
dimensional  crystal.  And  so,  computer  simulation  made  its  debut  in  sputtering 
research. 

A  very  important  early  computer  simulation  was  the  work  of  Gibson,  et  al 
[ref.  26]  in  1960.  They  developed  a  computer  model  which  simulated  the  motion  of 
primary  knock— on  atoms  in  a  copper  monocrystal  target.  Their  model  assigned  an 
arbitrary  kinetic  energy  and  direction  to  the  struck  atom  as  a  means  to  simulate  the 
collision  effect.  Atomic  interactions  were  treated  as  binary  collisions  and  the 
resulting  motion  obtained  through  Newtonian  mechanics.  Their  results 
demonstrated  the  ability  of  computer  simulations  to  isolate  elementary  sputtering 
prcesses  for  investigation.  Two  years  later  Robinson  and  Oen  [ref.  27],  utilizing  a 
similar  code,  discovered  the  channeling  effect.  What  makes  this  discovery  so 
important  is  the  fact  that  channeling  had  not  yet  been  observed  in  the  laboratory. 
For  the  first  time,  a  computer  simulation  had  predicted  an  important  dynamic 
sputtering  mechanism.  Channeling  was  later  verified  experimentally. 

In  1967,  Harrison,  Levy,  and  Effron  [ref.  28]  simulated  the  bombardment 
of  a  copper  monocrystal  by  argon  ions  using  repulsive  pair  potential  interactions. 
Their  findings  demonstrated  that  the  bulk  of  the  sputtering  yield  came  from  the 
first  three  layers  of  the  top  surface.  This  model  was  later  refined  to  include 


6 


attractive  potential  functions  for  the  target  [refs.  29—30).  This  potential  which 
includes  attractive  and  repulsive  interactions  produced  a  dynamically  stable  crystal. 
In  197S  Garrison.  Wi nograd  and  Harrison  [ref.  31]  published  the  result  of  a 
comprehensive  simulation  study  of  atomic  and  molecular  ejection  from  a  copper 
crystal  with  adsorbed  oxygen  atoms.  The  simulation  provided  a  detailed  description 
of  the  ejection  mechanism.  The  results  also  allowed  them  to  determine  whether 
molecules  were  ejected  as  clusters  or  whether  they  combined  together  after  leaving 
the  target's  surface.  They  also  determined  the  effects  of  adatom  placement  on 
molecule  formation. 

In  its  short  history  computer  simulation  has  successfully  predicted  and 
explained  various  specific  sputtering  events.  Results  of  simulations  have  been  both 
praised  and  strongly  criticized  by  some  theoreticians  who  argue  that  simulations 
neglect  important  factors.  Nevertheless,  the  value  of  computer  simulation  has  been 
recognized  bv  the  overall  scientific  communitv  and  it  is  now  considered  an 
important  research  area  in  the  field  of  sputtering. 

2.  General  Concepts  and  Development 

Most  re?l  systems  are  so  complicated  that  a  complete  analytical 
description  is  virtually  impossible  to  achieve,  unfortunately,  sputtering  falls  into 
this  category  of  complexity. 

The  modelling  of  a  real  collision  cascade  such  as  a  sputtering  event  is  an 
attempt  to  describe  a  complex  system  in  terms  of  idealized  and  highly  simplified 
descriptors.  Invariably,  many  factors  must  be  neglected  in  order  to  simplify  the 
model  the  selected  set  of  descriptors  must  adequately  determine  the  system's 
relevant  behavior.  Once  this  is  achieved,  the  result  is  a  simulation  of  a  real  system. 
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Computer  simulations  of  atomic  collisions  allow  researchers  to  concentrate 
on  the  basic  physics  of  a  system  without  the  analytic  constraints  of  statistical 
theories.  It  is  a  tool  that  can  be  used  to  test  the  applicability  of  a  theory. 

There  are  two  major  approaches  used  in  computer  simulations.  These  are 
the  time-step  and  discrete  event  models.  The  discrete  event  model  proceeds  from 
event  to  event  skipping  over  their  time  separation.  It  maintains  a  list  of  possible 
future  events  from  which  it  determines  the  next  event.  This  list  is  constantly  being 
updated.  This  model  works  well  when  the  events  are  sufficiently  separated  in  time. 
In  contrast,  the  time— step  model  advances  the  clock  for  a  short  interval  first,  and 
then  computes  all  the  interactions  that  occur  in  that  time,  updating  all  processes 
before  the  next  time  step.  This  model  is  the  obvious  choice  for  systems  with 
simultaneous  interactions.  The  time— step  logic's  chief  drawback  is  that  it  requires 
significantly  more  computer  running  time  than  the  discrete  event  program  logic. 

In  sputtering  applications,  these  two  approaches  lead  to  the  two  principal 
methods  of  simulation;  the  binary  collision  (BC)  and  the  multiple  interaction  (MI) 
program  logics. 

The  BC  method  [ref.  32]  uses  the  discrete  event  model.  The  basic 
assumption  of  this  type  of  simulation  is  that  each  particle  interacts  with  only  one 
other  particle  (usually  assumed  stationary)  at  a  time.  These  simulations  are 
inherentelly  restricted  to  be  linear  calculations. 

The  multiple  interaction  simulation  uses  the  time-step  model  and 
Newton's  equations  of  motion  are  numerically  solved  for  many  particles.  This  is  the 
method  used  in  this  thesis  investigation. 
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C.  SPUTTERING  OF  LIQUIDS 

The  sputtering  of  liquid  targets  has  gained  a  great  deal  of  attention  in  the 
scientific  community.  Experimental  results  of  sputtering  have  shown  that  solid 
targets  may  undergo  a  phase  transition  to  liquid  under  ion  bombardment. 
Therefore,  there  is  a  practical  as  well  as  scientific  interest  to  understand  and  be 
able  to  model  the  characteristics  of  a  liquid  surface. 

1.  Summary  of  Experimental  Results 

In  1962.  Whener  et  al  [ref.  33]  reported  their  measurements  on  tin 
sputtered  by  argon  ions.  Their  results  showed  a  40%  larger  yield  for  liquid  tin  at  an 
ion  energy  of  0.2  KeY  and  a  6%  smaller  yield  for  an  ion  energy  of  0.4  KeV.  This 
was  the  earliest  evidence  of  an  energy  dependent  yield  profile  with  a  maximum 
sputtering  yield  peak.  Eight  years  later  Krutenat  and  Panzera  [ref.  34]  conducted  a 
similar  sputtering  experiment  on  solid  and  liquid  tin.  Their  results  were  in  general 
agreement  with  those  of  Whener,  showing  that  the  solid  yields  below  0.2  KeV  is 
twice  that  of  the  liquid  and  that  a  crossover  occurs  at  0.375  KeV  above  which  the 
liquid  yield  is  only  15%  higher.  Their  results  also  implicated  surface  dynamics  in 
sputtering  mechanisms. 

In  1982.  Cooper  and  Hurst  [ref.  35]  bombarded  liquid  and  solid  indium 
with  argon  ions  and  discovered  that  above  0.107  KeV  the  liquid  yield  was  10% 
higher  than  the  solid.  They  noted  that  the  shape  of  the  yield  versus  ion  energy  for 
both  the  liquid  and  solid  were  similar.  They  attributed  the  higher  yield  of  the 
liquid  to  its  lower  binding  energy. 

Dumke  et  al  [ref.  36]  published  their  results  on  the  sputtering  of  a 
gallium— indium  eutectic  alloy  (surface  monolayer  greater  than  94%  indium)  in  the 
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liquid  phase.  Their  findings  indicated  a  ratio  indium  to  gallium  sputtering  yield 
which  was  28  times  greater  than  expected  from  the  stoichiometry.  They  were  able 
to  conclude  that  85%  of  the  sputtered  atoms  came  from  the  surface  monolayer. 

2.  Historical  Summary  of  Liquid  Models 

In  order  to  study  the  effects  of  sputtering  on  a  liquid  surface  using 
computer  simulations,  one  first  needs  to  develop  a  reasonable  model  of  a  liquid. 

One  very  early  model  of  a  simple  liquid  was  developed  by  Bernal  in  1960 
[ref.  37],  He  stipulated  that  the  structure  of  a  simple  liquid  can  be  determined 
simply  by  volume  exclusion.  He  then  proposed  a  zeroth  order  model  in  which  the 
atoms  were  considered  hard  spheres  and  their  local  structure  determined  by  the 
restriction  that  no  two  atoms  can  approach  each  other  by  less  than  one  atom 
diameter.  Early  Bernal  models  were  constructed  [refs.  38,39]  by  pouring  thousands 
of  steel  balls  into  a  deformable  container  (originally  a  football  bladder).  The 
ensemble  was  then  bound  with  rubber  strips  and  kneaded  to  maximize  its  density. 
The  contents  were  then  fixed  by  pouring  melted  wax  into  the  container  through 
holes  and  then  allowed  to  solidify.  The  coordinates  of  each  sphere  were  then 
individually  and  painstakingly  measured  (probably  by  one  of  his  postgraduate 
students).  This  seemingly  ad  hoc  procedure  yielded  a  maximum  hard— sphere 
non— crystalline  packing  of  0.6366.  Amazingly,  more  recent  attemts  to  reproduce 
this  structure  using  modern  computational  techniques  have  at  best  produced 
indistinguishable  results. 

In  1968,  Verlet  [ref.  40]  published  the  results  of  his  computer  experiments 
on  classical  fluids.  His  was  one  of  the  earliest  attempts  at  modelling  the  liquid 
structure  through  a  computer  simulation.  His  model  used  a  Lennard— Jones 
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potential  as  the  means  of  particle  interaction  from  which  the  equations  of  motion 
were  solved.  His  study  offered  evidence  that  correlation  functions  at  high  density  is 
due  to  the  geometrical  effects  of  a  strong  repulsion  in  the  potential.  He  showed  that 
the  same  behavior  can  be  obtained  by  the  approximate  solution  of  the  hard  sphere 
problem  and  that  the  diameter  of  the  hard  spheres  is  the  only  parameter  of  the 
theory  (in  agreement  with  the  simplistic  Bernal  model). 

In  1982,  Miranda  and  Torra  [ref.  41]  obtained  good  results  for  the  liquid 
structure  and  for  the  self— diffusion  constant  with  a  computer  simulation  using 
functions  derived  from  a  local  pseudopotential  and  various  other  dielectric  functions. 

3.  Computer  Simulations  Summary  of  Liquid  Sputtering 

In  1985,  Lo  et  al  [ref.  42]  generated  two  liquid  targets  consisting  of  603 
copper  atoms  and  then  subjected  them  to  an  argon  ion  beam.  The  two  targets 
differed  only  in  the  imposed  boundary  conditions.  One  target  used  a  box  boundary 
condition  requiring  particles  to  experience  pure  reflection  at  the  boundaries  and  the 
other  used  a  semi— periodic  boundary  condition  requiring  position  and  momentum 
periodicity  in  the  two  dimensions  defined  by  the  surface.  They  found  that  the 
semi— periodic  boundary  condition  results  were  in  better  agreement  with 
experimental  sputtering  results  and  concluded  that  such  a  target  better  represented 
the  free  surface  of  a  real  liquid.  The  target  was  warmed  to  liquid  temperatures  by 
assigning  velocities  to  each  atom  with  a  random  number  generator.  The  target  was 
then  allowed  to  equilibrate  for  a  few  picoseconds.  This  energy  imput  scheme  is 
particularly  important  to  this  research  effort  because  it  is  used  in  generating  one  of 
our  targets. 
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In  19S7.  Lo  et  al  [ref.  43]  again  conducted  a  computer  simulation  study  of 
collision  cascades  in  liquid  indium.  Their  results  suggested  that  the  detailed 
structure  of  the  target  surface  layer  is  very  important  in  the  sputtering  proccess. 

The  most  recent  simulation  of  liquid  sputtering  was  done  by  Morgan 
[ref.  44],  He  conducted  a  study  of  self-sputtering  using  a  stratified  liquid  metal 
surface  as  a  model.  He  observed  an  enhanced  low  energy  yield  which  fell  below 
those  of  other  models  for  higher  energies.  These  results  appear  to  be  in  general 
agreement  with  various  published  measurements  of  liquid  sputtering  yields 
[refs.  34-36], 
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II.  OBJECTIVES 


This  research  and  computer  simulation  effort  has  three  specific  objectives.  The 
first  objective  is  to  develop  a  reasonable  model  of  copper  in  the  liquid  phase  by 
using  a  heavily  modified  version  of  QDYN  (Quick  Dynamics),  a  molecular  dynamics 
computer  program  developed  by  Harrison  [ref.  45].  The  general  intent,  is  to  warm 
an  fee  (copper)  crystal  by  giving  each  atom  in  the  crystal  sufficient  energy  to  reach 
liquid  temperatures.  The  energy  given  to  each  atom  is  determined  stochastically. 
The  warmed  ensemble  is  then  allowed  to  interact  through  pair-wise  forces  until 
thermal  equilibrium  is  reached  at  the  desired  liquid  temperature.  The  resulting 
liquid  targets  will  be  compared  to  various  known  liquid— like  characteristics  such  as 
the  radial  distribution  function  for  copper  in  the  liquid  state  and  the  expected 
Maxwell— Boltzmann  (velocity  distribution)  function. 

Eight  different  liquid  targets  will  be  generated  using  two  different  methods  of 
warming  (adding  energy  to  the  atoms)  in  combination  with  four  different  variations 
of  dynamic  integration.  These  variations  are  the  permuted  combinations  of 
simulation  runs  with  or  without  a  reflective  boundary  condition  and,  or,  with  or 
without  the  requirement  to  update  the  nearest  neighbor  list  in  each  time  step. 

The  second  objective  is  to  compare  the  dynamic  behavior  of  the  two  different 
warming  methods.  One  of  these  warming  methods  is  a  variation  of  the  method  used 
by  Lo  et  al  [ref.  42],  where  the  energy  is  added  to  each  atom  by  assigning  each 
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a  random  velocity.  The  second  warming  method  was  one  suggested  by  the  late 
professor  Harrison  of  the  Naval  Postgraduate  School.  This  involves  warming  a 
target  by  randomly  displacing  each  atom  from  their  stable  positions  and  then 
allowing  then  to  interact  through  pair-wise  forces  until  equilibration  is  reached. 
Clearly,  each  method  starts  the  molecular  dynamics  with  a  different  form  of  excess 
energy.  One  method  starts  with  excess  kinetic  energy  and  the  other  starts  with 
excess  potential  energy.  Therefore,  the  ability  of  QDYN  to  equilibrate  the  targets 
from  such  varying  initial  energy  conditions  will  serve  as  further  evidence  for  the 
applicability  of  pair— potentials  in  molecular  dynamics. 

The  third  objective  is  to  use  the  resulting  liquid  targets  in  a  computer 
sputtering  experiment.  The  goal  is  to  obtain  sputtering  yields  which  are  in  general 
agreement  with  published  experimental  data  [refs.  33—36]  as  well  as  other  liquid 
sputtering  simulation  results  [refs.  42 — 44] . 
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III.  THEORY 


A.  LIQUIDS 

Matter  manifests  itself  in  three  states.  One  of  these  is  the  liquid  state  and  is 
the  least  understood  of  all  three.  Understanding  the  unique  descriptive 
characteristics  of  a  liquid,  especially  those  that  may  affect  the  interatomic  collision 
process,  such  as  density,  radial  distribution,  velocity  distribution  and  dynamic 
behavior,  is  of  great  importance  to  this  research  effort. 

1.  The  Density  of  Liquids 

The  essential  difference  between  a  solid  and  a  liquid  is  that  a  liquid  is  a 
much  more  disordered  state  [ref.  47].  A  disordered  ensemble  of  atoms  requires  more 
space  than  does  a  typical  closed-packed  and  well-ordered  crystal  structure. 
Consequently,  the  density  of  a  liquid  structure  can  be  expected  to  have  lower 
density  than  that  of  a  solid.  In  the  case  of  copper,  the  density  has  been 
experimentally  measured  for  a  wide  range  of  temperatures  in  the  liquid  phase  and  is 
found  to  linearly  decrease  with  increasing  temperatures.  Cahill  and  Kisherbaum 
[ref.  4S]  have  obtained  experimental  density  measurements  for  copper  from  the 
melting  point  (1356  K)  to  2500  K.  Their  results  indicate  that  a  plot  of  the  density 
of  copper  versus  temperature  in  the  liquid  phase  can  be  accurately  described  by  a 
straight  line  of  negative  slope  given  by, 

pl  =  9.077  -  (8.006xl0'4  TJ 

where  Tj  is  the  liquid  temperature  of  copper  in  degrees  Kelvin  and  pl  is  the  density 
of  copper  in  grams  per  cubic  centimeter. 
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2.  The  Radial  Distribution  of  Liquids 

The  liquid  state  imposes  spatial  restrictions  on  the  possible  ways  that 
atoms  can  arrange  to  form  amorphous  structures.  The  nearest  neighbors  to  a  given 
atom  cannot  move  too  far  away  because  of  the  repulsion  force  exerted  by  the  other 
particles  in  the  system.  Due  to  these  restrictions  the  liquid  state  has,  on  average,  a 
very  distinctive  radial  distribution  as  described  by  Azaroff  [ref.  49].  This  short 
range  order  of  liquid  makes  it  possible  to  describe  atomic  structure  in  terms  of  the 
average  density  of  atoms  per  unit  volume  pQ  and  the  actual  density  p{ r)  a  distance  r 
away  from  an  atom  placed  at  the  origin. 

In  monatomic  liquids  such  as  copper,  the  radial  density  is  spherically 
symmetric  and  can  be  described  as  the  number  of  atoms  per  unit  volume  contained 
in  a  spherical  shell  of  radius  r  and  of  thickness  dr  around  one  atom  designated  as  the 
origin.  The  volume  of  such  a  shell  is  the  difference  between  the  volumes  of  two 
spheres  whose  radii  are  r  and  r+dr.  The  density  of  atoms  in  such  a  shell  is  given  by 
4  7rr2/?(  r  )dr. 

The  radial  distribution  function  for  the  atoms  is  defined  as. 

G(r)  -  47rr2p(r). 

The  radial  distribution,  however,  cannot  be  experimentally  measured. 
Instead,  the  structure  factor  i(s).  defined  below,  is  obtained  through  x— ray  or 
neutron  diffraction  experiments  from  which  the  radial  distribution  is  calculated. 

The  radial  distribution  function  G(r)  is  related  to  the  structure  factor  i(s) 
by  the  following  Fourier  transformation  [ref.  50]: 

G(r)  =  47rr2p(r)  =  4;rr 2p  +  —  f  s  i(s)  sin(rs)  ds 

0  7T  J 
O 

where  pQ  is  the  average  density  of  the  ensemble. 
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The  problem  with  this  transformation  is  that  the  values  for  i(s)  must  be 
known  through  oc  .  To  avoid  this  problem,  the  method  of  Zei  and  Steffan  [ref.  51] 
can  be  used  to  extend  the  range  of  measured  i(s)  without  manipulating  the 
experimental  values.  This  is  the  method  used  by  Eder  et  al  [ref.  46]  to  obtain  the 
radial  distribution  of  liquid  copper  at  1393K  and  1833K.  This  is  the  radial 
distribution  against  which  the  results  of  this  simulation  will  be  compared. 

In  computer  simulations,  the  radial  distribution  is  very  easily  obtained  by 
tallying  the  actual  separations  of  each  atom  from  every  other  atom  in  the  ensemble. 
This  is  the  method  used  in  our  simulation.  The  distribution  is  normalized  by 
dividing  each  bin  total  by  the  area  of  a  shell  whose  radius  is  equal  to  the  minimum 
distance  of  the  bin.  The  reason  for  doing  this  is  to  better  compare  the  distribution 
of  a  very  small  ensemble  with  that  of  a  real  system  whose  size  is  infinitely  larger  by 
comparison. 

3.  The  Velocity  Distribution  of  Liquids 

The  general  collision  behavior  of  liquid  particles  as  they  interact  through 
pair-wise  forces  is  not  very  different  from  that  of  a  gas.  Therefore,  we  can  expect  a 
velocity  distribution  of  a  liquid  to  be  similar  to  the  velocity  distribution  of  a  gas. 

The  distribution  of  molecular  speeds  in  a  large  sample  of  gas  varies  over  a 
wide  range  of  magnitude  and  has  a  characteristic  distribution  which  depends  on 
temperature.  The  first  expression  for  the  velocity  distribution  of  a  gas  was  derived 
out  by  Clerk  Maxwell.  According  to  Maxwell,  a  sample  of  gas  containing  N 
molecules  has  a  distribution  of  velocities  given  by  [ref.  52]: 


N(v)  =  4trNt  m/2trkT  2  v2  e  (mv'2/2kT) 


where  N(v)dv  is  the  number  of  molecules  having  speeds  between  v  and  v+dv,  T  is 
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the  absolute  temperature,  k  is  the  Boltzmann's  constant,  Nt  is  the  total  number  of 
molecules  and  m  is  the  mass  of  the  molecule. 

4.  Equilibrium  and  the  Equipartition  of  Energy 

The  equipartition  of  energy  theorem  states  that  when  the  number  of 
particles  is  large  and  Newtonian  mechanics  hold,  all  the  energy  terms  have  the  same 
average  value,  and  that  the  average  value  depends  only  on  the  temperature  [ref.  52]. 
This  means  that  the  available  energy  distributes  itself  in  equal  shares  to  each  of  the 
independent  ways  in  which  molecules  can  absorb  energy. 

The  energy  given  to  the  atoms  by  the  warming  scheme  in  the  simulation 
will  distribute  itself  into  half  kinetic  energy  and  half  potential  energy.  The  energy 
going  into  kinetic  will  further  distribute  itself  equally  among  the  x,  y  and  z 
components  of  kinetic  energy.  The  energy  is  distributed  through  the  many  collisions 
that  the  atoms  experience  as  the  ensemble  comes  to  equilibrium.  The  target  is  at 
equilibrium  when  the  total  added  energy  is  completely  partitioned. 

A  plot  of  the  total  kinetic  energy  versus  time  should  resemble  a  damped 
sinusoid,  oscillating  about  a  value  equal  to  half  of  the  total  added  energy.  A  plot  of 
the  x,v  and  z  components  of  kinetic  versus  time  should  show  the  components 
becoming  equal  as  time  progresses  and  their  value  approaches  one  sixth  of  the  total 
added  energy. 

B.  SPUTTERING 

1.  General  Concepts 

Sputtering  is  the  erosion  of  a  target  surface  through  the  removal  of  atoms 
by  the  action  of  incident  energetic  particles.  Sputtering  occurs  in  nature  when  a  hot 
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plasma,  such  as  that  found  in  a  lightening  bolt,  comes  in  contact  with  a  solid 
surface.  Sputtering  can  be  produced  in  the  laboratory  by  subjecting  a  target  surface 
to  a  plasma  or  a  particle  ion  beam. 

The  most  obvious  measure  of  erosion  effects  is  through  the  sputtering 
yield,  Y,  defined  as  the  average  number  of  atoms  removed  per  incident  particle. 
The  yield  is  a  function  of  target  structure,  target  atom  mass,  incident  beam 
alignment  (relative  to  target  structure),  energy  of  incident  particles  and  the  mass  of 
the  incident  particle.  Sputtered  particles  can  leave  the  target's  surface  with  a  broad 
distribution  of  energy,  charge  state  and  exit  angles.  In  order  to  obtain  reproducible 
experiments,  the  following  conditions  must  be  met  [ref.  53]: 

•  The  target  surface  must  be  clean. 

•  The  gas  pressure  must  be  low  enough  such  that  the  mean  free  path  of  ions 
and  sputtered  atoms  is  large. 

•  The  ion  current  density  must  be  high  and  the  background  pressure  low  so 
that  formation  of  surface  layers  is  prevented  during  the  experiment. 

•  The  ions  must  strike  the  target  at  a  known  angle. 

•  The  energy  spread  of  the  incident  beam  must  be  small. 

•  The  ionizing  conditions  in  the  ion  source  should  minimize  the  production  of 
multiple  charged  species;  the  atoms  must  be  uniformly  charged  and 
separated. 

•  The  lattice  orientation  of  monocrystalline  targets  or  texture 
of  polvcrystal  must  be  known. 

Sputtering  theory  has  four  major  variants;  collisional,  thermal,  electronic, 
and  exfoliational.  In  this  research  effort,  we  are  primarily  interested  in  the 
collisional  behavior  of  sputtering. 
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2.  Collisional  Description 


The  collisional  description  of  sputtering  views  the  bombarding  ions  as 
energetic  projectiles  which  deposit  kinetic  energy  and  momentum  upon  striking  a 
target  atom.  This  action  is  analogous  to  an  atomic  size  game  of  billiards. 

There  are  two  main  collisional  processes;  prompt  collisional  and  slow 
collisional.  Their  name  arises  from  their  respective  collision-event  time  duration. 
Once  again  we  limit  the  scope  of  the  research  by  concentrating  on  the  prompt 
collisional  process. 

The  prompt  collisional  process  has  a  duration  time  of  approximately  500 
femptoseconds  following  impact.  The  general  process  occurs  in  the  following 
manner;  an  incoming  particle  collides  with  a  target  atom  and  transfers  some  of  its 
kinetic  energy  to  it.  If  the  energy  received  by  the  target  atom  exceeds  the  binding 
energy  of  its  lattice  site,  a  primary  knock— on  atom  is  created.  This  knock-on  atom 
can  travel  through  the  target's  lattice  colliding  with  other  atoms  and  transferring 
some  of  its  kinetic  energy  and  momentum.  A  near  surface  atom  is  sputtered  if  its 
kinetic  energy  has  a  component  normal  to  the  surface  which  is  larger  than  the 
surface  potential  energy  barrier. 

There  are  four  basic  collisional  sputtering  descriptions;  rebound,  recoil, 
reflection,  and  direct  or  reflected  recoil,  these  are  summarized  in  turn  in  the 
following  sub— sections. 

a.  Rebound  Sputtering 

In  the  rebound  sputtering  sequence,  an  incident  ion  strikes  an  atom 
in  the  first  atom  layer.  The  atom  initially  receives  a  component  of  kinetic  energy 
normal  to  the  target  surface  and  into  the  target.  The  atom  then  collides  and  recoils 
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off  an  atom  in  the  second  layer  reversing  the  surface  normal  component  of  kinetic 
energy.  The  atom  then  escapes  the  surface  through  the  first  atom  layer  as 
illustrated  by  Figure  1  of  Appendix  A. 

This  mechanism  applies  only  to  adatoms  which  are  much  lighter 
than  atoms  of  the  substrate  and  assumes  that  atoms  recoil  without  energy  loss. 

b.  Recoil  Sputtering 

In  the  recoil  sputtering  process,  the  first  atom  struck  by  the 
bombarding  particle  is  not  the  ejected  atom.  The  bombarding  particle  imparts 
kinetic  energy  and  momentum  to  an  atom  in  the  first  atom  layer  in  a  collision 
process  similar  to  the  one  described  in  rebound  sputtering.  The  atom  then 
undergoes  a  sequence  of  collisions  with  other  target  atoms  as  it  plows  inwards.  The 
atom's  initial  inward  component  of  kinetic  energy  is  eventually  turned  back  towards 
the  surface  through  multiple  recoil  collisions.  The  atom  finally  strikes  another  atom 
from  below  causing  its  ejection  from  the  surface.  Often,  only  two  collisions  are 
required  as  illustrated  by  Figure  2  of  Appendix  A. 

c.  Reflection  Sputtering 

In  reflection  sputtering,  the  incoming  ion  is  backscattered  by  a 
target  atom  below  layer  one.  The  reflected  ion  strikes  an  atom  in  layer  one  from 
below  causing  its  ejection  (see  Figure  3  in  Appendix  A).  This  process  is  extremely 
rare  due  to  the  space  and  angle  restrictions  imposed  by  the  lattice  geometry. 

d.  Sputtering  by  Direct  or  Deflected  Recoil 

In  direct  or  deflected  recoil  sputtering,  an  atom  in  any  near-surface 
layer  is  struck  by  the  incoming  particle  at  grazing  incidence.  The  struck  atom  is 
deflected  by  neighboring  atoms,  causing  its  ejection  from  the  target  surface.  This 
process  is  illustrated  by  Figure  4  of  Appendix  A. 
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The  four  collisional  sputtering  mechanisms  described  above  are 
simplified  and  highly  idealized.  What  most  often  occurs  in  a  sputtering  sequence, 
however,  is  the  formation  of  collisional  cascades.  When  the  first  atom  in  a  sequence 
is  struck  by  the  bombarding  particle,  it  will  undergo  a  series  of  collisions  throughout 
the  crystal  structure  until  the  atom  exhausts  its  kinetic  energy  or  it  is  ejected  from 
the  target  structure.  These  multiple  collision  sequences  allow  many  paths  through 
which  the  initial  ion  energy  can  reach  the  surface  of  the  target  and  cause  sputtering 
or  reach  deep  target  atom  layers  and  dissipate  its  energy  as  heat. 

As  a  collision  cascade  develops,  it  may  cause  the  sputtering  of  atoms 
through  any  of  the  mechanisms  described  above  or  a  combination  thereof.  However, 
it  is  also  possible  for  a  cascade  to  produce  no  ejections  at  all.  Such  a  cascade  is 
illustrated  by  Figure  5  of  Appendix  A. 

C.  THE  COMPUTER  MODELS 
1.  Experimental  Approach 
a.  General  Description 

The  general  approach  to  this  study  is  to  first  obtain  several 
reasonable  liquid  targets  and  then  use  the  the  best  resulting  target  in  a  sputtering 
run.  The  first  step  is  to  generate  a  fcc(lll)  crystal  lattice.  This  is  done  with  the 
subroutine  Fill  (ref.  54].  Next,  the  atom  ensemble  is  energized  using  a  stochastic 
scheme.  In  this  simulation,  there  are  two  distinct  methods  used  to  randomly  assign 
energy  to  individual  atoms.  These  will  be  individually  discussed  in  later  sections. 

After  the  excess  energy  is  given  to  the  target,  the  dynamic  part  of 
the  liquid  simulation  starts,  and  the  atoms  are  allowed  to  interact  through 
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pair-wise  forces  until  they  reach  thermal  equilibrium.  This  part  of  the  simulation  is 
done  with  a  modified  version  of  QDYN  [ref.  45].  In  sputtering  simulations  using 
QDYN,  atoms  are  only  assumed  to  move  after  collision  with  a  moving  atom,  in 
generating  the  liquid  all  the  atoms  are  designated  as  moving  atoms  from  time-step 
one. 

b.  Assumptions 

The  most  significant  assumption  is  that  the  motion  of  particles  in  a 
liquid  state  can  be  approximated  by  pair— potential  derived  forces  and  that  the 
particles  obey  Newtonian  Mechanics.  Another  assumption  (for  sputtering)  is  that 
the  atoms  in  a  liquid  target  move  so  slow,  compared  to  the  ion  and  collision  cascade 
atoms,  that  they  can  be  considered  static.  The  average  velocity  of  liquid  atoms  at 
1700K  is  about  670  m/sec.  The  average  collision  cascade  atom  is  at  least  10  times 
faster  and  the  ion  is  almost  3000  times  faster.  A  one  KeV  Argon  ion  moves  at  3030 
km/sec.  So  during  a  typical  cascade  lasting  300  fs  the  average  liquid  atom  moves 
2xl0'2  m,  approximately  one  half  of  a  lattice  unit. 

2.  Liquid  Target  Generation 

The  basic  target  consists  of  1445  atoms  originally  placed  at  the  lattice 
positions  of  an  fee  crystal  with  a  111  plane  orientation.  The  atom  positions  are 
defined  in  terms  of  lattice  units.  The  lattice  unit  for  copper  is  consistently  used 
throughout  the  simulation  as  the  basic  unit  of  length.  This  unit,  as  well  as  other 
copper  constants  used  in  the  simulation,  can  be  found  in  Table  1  of  Appendix  B. 
The  dimensions  of  the  crystal  are  19  x  8  x  17  lattice  units. 
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The  target  is  aligned  in  such  a  way  that  an  incoming  ion  would  strike  the 
target  in  the  111  plane  and  travel  in  the  positive  y  direction  as  illustrated  by 
Figure  6  of  Appendix  A. 

The  atoms  are  warmed  using  two  techniques  as  explained  in  the  following 
sub— sections. 

a.  Warming  by  Impulse 

This  method  of  warming  atoms  is  a  variant  of  the  method  used 
by  Lo  et  al  [ref.  42].  The  general  idea  is  to  give  each  atom  an  initial  impulse  whose 
direction  and  n  ignitude  is  determined  by  a  random  number  generator. 

The  warming  is  accomplished  in  two  steps  and  each  uses  a 
different  random  number  generator.  The  magnitude  of  velocity  for  each  atom  is 
obtained  from  the  product  of  a  multiplicative  constant  (the  speed  of  interest)  and  a 
set  of  random  numbers  having  a  Gaussian  distribution  about  unity  with  a  standard 
deviation  of  one.  The  random  numbers  were  scaled  to  go  from  zero  to  two.  This 
was  done  to  avoid  large  velocities.  The  multiplicative  velocity  constant  is  chosen  to 
be  the  velocity  corresponding  to  twice  the  desired  kinetic  energy.  This  is  done  to 
allow  for  the  equipartition  of  half  of  the  added  energy  into  potential  energy 
elements.  In  our  case,  the  desired  temperature  is  1462  K.  Therefore,  the  velocity 
corresponding  to  twice  the  temperature  is  1067.7  meters  per  second.  This  is  the 
multiplicative  constant  used  in  the  generation  of  random  speeds.  The  individual 
atom  velocities  are  the  products  of  this  speed  with  a  series  of  random  numbers. 

The  direction  of  each  atom  must  be  uniformly  distributed  in  3— D 
space.  The  velocity  directions  are  obtained  from  the  sines  and  cosines  of  a  spherical 
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coordinate  system.  The  sines  and  cosines  are  obtained  from  a  random  number 
generator  which  produces  uniform  random  numbers  between  zero  and  one. 

The  actual  warming  is  done  by  a  subroutine  called  "warm",  excerpt 
of  which  is  shown  below, 


SUBROUTINE  WARM 
PARAMETER(LTMX=1854) 

DIMENSION  RAN(LTMX),VRAN(LTMX*5) 


DO  300  1=1, LT 
R1=VRAN(I) 

R2=VRAN(I+(2*LT)) 

R3=VRAN(I+(3*LT)j 

R4=VRAN(I+(4*LT)) 

R5= VRAN  (1+  (5*LT)) 

VELMAG=VMAG*RAN(I) 

SIGN1=+1.0D0 

IF(R5.LT.0.5D0)SIGN1=— 1.0D0 
SIGN2=+1.0D0 

IF(R5.LT.0.5D0)SIGN2=— I.0D0 
COSPHI=SIGN2*Rl 

SINPHI=DSQRT(  1 .0D0-COSPHI*COSPHI) 
SINTHE=SIGN1*(2.0D0*R3*R4)/(R3*R3+R4*R4) 
COSTHE=(R3*R3— R4*R4)/(R3*R3+R4*R4) 
VX(I)=VELMAG*(SINPHI*COSTHE) 
VY(I)=VELMAG*(SINPHI*SINTHE) 
VZ(I)=VELMAG*COSPHI 
300  CONTINUE 


CONTINUE 

RETURN 

This  do— loop  goes  from  1  to  LT,  where  LT  is  the  total  number  of 
atoms.  VELMAG  is  the  velocity  magnitude  of  each  atom  I,  obtained  from  the 
product  of  a  velocity  multiplicative  constant,  VMAG,  and  a  random  number 
RAN(I).  The  vector  RAN(I)  is  a  set  of  random  numbers  with  a  gaussian 
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distribution  about  unity  and  with  a  standard  deviation  of  one.  These  numbers  are 
not  allowed  to  be  greater  than  two  or  less  than  zero  to  avoid  large  velocities. 
R1-R5  are  uniformly  distributed  sets  of  random  numbers  obtained  from  the  vector 
VRAN(I).  R2  and  R5  are  used  to  define  the  sign(+/— )  for  SINTHE  and  COSPHI. 
COSPHI,  SINPHI.  SINTHE,  and  COSTHE  are  the  sines  and  cosines  of  the  usual 
spherical  coordinate  directional  angles  y  and  a.  VX(I),  VY(I),  and  VZ(I)  are  the 
resulting  components  of  velocity  for  atom  I. 

The  liquid  simulation  program  using  this  type  of  warming  technique 

is  called  QLV. 

b.  Warming  by  Displacement 

This  warming  scheme  gives  added  potential  energy  to  the  ensemble 
by  individually  displacing  each  atom  away  from  their  positions  of  minimum 
potential.  The  displacement  is  assigned  to  each  atom  by  adding  multiples  of  the 
thermal  amplitude  to  the  x,v  and  z  components  of  position.  The  individual 
component  displacement  is  the  product  of  the  thermal  amplitude  and  a  random 
number.  The  random  number  generator  used  produces  numbers  which  are  Gaussian 
distributed  about  zero  and  with  a  standard  deviation  of  one. 

The  thermal  amplitude,  TaD)p,  is  obtained  from  the  relation, 


where  TT^  is  the  mean  square  of  total  displacements  of  an  atom  from  the  average 
position.  This  is  a  value  which  can  be  experimentally  obtained. 

The  existing  experimental  measurements  of  u-2^  are  for  temperatures 
well  bebw  the  liquid  phase.  For  this  reason,  the  warming  must  be  done  in  cycles. 


26 


Each  atom  is  displaced  several  times  until  a  suitable  liquid  temperature  is  achieved. 
The  first  warming  cycle  consists  of  displacing  all  atoms  from  their  equilibrium 
positions  by  adding  to  each  component  of  position  a  random  multiple  of  a  known 
thermal  amplitude.  The  following  warming  cycles  further  displaces  each  atom  from 
their  last  position  using  the  same  principle.  The  warming  behavior  of  such 
sequential  cycling  is  that  a  maximum  average  temperature  is  reached  after  several 
cycles.  Further  displacements  cause  the  total  energy  to  oscillate  about  some 
average. 

The  thermal  amplitude  corresponding  to  300K  was  chosen  because  it 
was  an  actual  data  point  in  the  experimental  results  of  Singh  and  Sharma  [ref.  55]. 
The  corresponding  thermal  amplitude  used  in  the  simulation  is  0.0782  lattice  units. 
The  liquid  simulation  program  using  this  type  of  warming  technique  is  called  QLP. 
c.  Achieving  Thermal  Equilibrium 

Thermal  equilibrium  is  achieved  by  allowing  the  atoms  to  interact 
with  each  other  through  pair— wise  forces.  The  target  is  considered  to  be  in  thermal 
equilibrium  when  the  plot  of  kinetic  energy  versus  time  settles  at  half  of  the  total 
added  energy  and  the  x,y  and  z  components  settle  at  equal  values.  Since  the  added 
energy  in  the  impulse— warmed  program  starts  with  all  kinetic  energy  and  the 
displacement— warmed  program  starts  with  all  potential  energy,  the  plots  of  kinetic 
energy  versus  time  for  QLP  and  QLV  resemble  mirror  images.  These  plots  are 
similar  to  damped  sinusoids  which  approach  the  value  of  half  of  the  total  added 
energy  as  the  target  comes  to  equilibrium. 
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d.  Boundary  Conditions 

Simulations  were  run  with  and  without  boundary  conditions.  The 
boundary  conditions  requires  the  atoms  to  experience  pure  reflection  at  the 
boundaries.  The  target  boundaries  for  reflection  are  the  uniformly  expanded  solid 
crystal  dimensions  corresponding  to  a  liquid  density.  The  density  of  copper  is 
obtained  from  the  experimental  measurements  of  Cahill  and  Kirshenbaum  [ref.  48]. 

The  simulation  programs  were  also  run  without  boundary  conditions 
to  investigate  the  extent  of  the  overall  system  expansion  due  to  the  interatomic 
interactions. 

3.  Sputtering  Program 

The  computer  program  used  for  sputtering  is  SDYN88A,  and  it  is  a 
variant  of  QDYN  developed  by  Harrison  [ref.  45].  The  program  uses 
multiple -interaction  (MI)  logic  in  a  time-step  approach. 

The  following  actions  occur  in  each  time— step; 

•  Summation  of  pair-wise  forces  for  each  atom. 

•  Calculation  of  new  velocities  and  positions. 

•  Movement  of  atoms  to  their  new  positions. 

•  Test  energy  conservation. 

The  forces  affecting  individual  atoms  are  not  computed  until  a  collision 
occurs.  This  initial  collision  turns  the  atoms  on  (designates  them  as  moving  atoms) 
and  puts  them  on  a  moving-atom  list.  Once  in  this  list,  the  atom's  interactions  are 
computed  each  time  step. 

The  moving  atoms  maintain  a  list  of  the  atoms  with  which  they  can 
interact.  This  list  is  called  the  nearest  neighbor  list.  The  process  of  updating  this 
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list  requires  considerable  computing  time.  For  this  reason,  it  is  updated  every  five 
time  steps. 

The  forces  are  computed  by  solving  Newton's  equations  of  motion.  A 
predictor— corrector  integration  scheme  adjusts  the  time  increment  of  each  time  step 
based  on  the  fastest  atom.  The  integration  technique  is  unique  in  that  position  and 
velocity  are  calculated  for  the  same  instant  of  time. 

The  program  w-as  originally  designed  to  study  the  effects  of  sputtering  in 
solid  targets.  Consequently,  the  original  logic  did  not  calculate  the  new  velocities 
and  positions  of  every  atom  each  time— step.  Instead,  only  the  moving  atoms  were 
updated.  In  our  version,  however,  all  of  the  target  atoms  have  velocity  from  the 
start  of  dynamic  integration.  The  initial  atom  positions  and  velocities  are  read 
from  an  input  file  created  by  the  liquid  simulation  programs  QLP  and  QLV. 

In  order  to  study  the  validity  of  using  a  static  target  in  liquid  sputtering 
simulations,  we  have  allowed  for  the  simulation  to  proceed  with  either  the  atoms 
turned  "on"  or  "off".  When  the  atoms  are  turned  on,  the  simulation  assumes  every 
atom  in  the  target  is  moving.  When  the  atoms  are  turned  off,  the  individual  atom 
velocities  are  zeroed  until  they  undergo  a  collision  with  a  fast  atom  or  ion.  Such 
static  targets  are  the  equivalent  of  an  amorphous  solid  with  liquid  densities.  The 
computing  time  for  a  run  with  the  atoms  turned  on  is  about  five  times  greater.  The 
justification  for  using  a  static  target  is  that  the  velocity  of  a  IKeV  ion  is  much 
faster  than  that  of  an  atom  in  the  liquid  state  (about  3000  times  faster  for  1700  K). 
The  average  velocity  of  a  Cu  atom  in  the  liquid  state  is  about  800  m/sec.  The 
collision  cascade  atoms  are  about  10  times  faster  than  liquid  atoms. 
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The  ion  impact  points  are  determined  by  a  grid  of  300  irregularly  spaced 
points.  The  grid  is  one  lattice  unit  squared.  T he  grid  is  overlaid  near  the  center  of 
the  target.  For  our  target  the  grid  reference  point  was  chosen  as  x=6.25  LU  and 
y=  11.5  LU.  The  impact  point  coordinates  are  read  in  from  the  input  File  and  added 
to  the  reference  position. 

D.  POTENTIAL  FUNCTIONS 
1.  Background 

If  we  take  two  atoms  originally  separated  by  a  large  distance  and  slowly 
bring  them  together,  the  atoms  will  eventually  start  to  experience  an  attractive 
force.  This  attraction  can  be  simulated  with  an  attractive  potential  function  that 
has  the  characteristics  of  a  well.  As  the  particles  get  closer,  the  potential  decreases 
until  the  point  where  the  force  changes  from  attractive  to  repulsive.  This  point 
corresponds  to  the  equilibrium  separation.  As  the  pair  gets  closer  together,  the 
repulsive  force  quickly  increases.  This  is  known  as  the  hard  collision  range.  The 
potential  in  this  range  has  the  characteristics  of  a  wall. 

In  the  simulation,  the  molecular  dynamic  interactions  are  caused  by- 
two— body  forces.  The  basic  assumption  is  that  the  forces  and  potential  energies 
depend  solely  on  the  physical  properties  of  the  interacting  particles  and  their 
internuclear  separation.  The  individual  pair-wise  interactions  that  each  atom 
experiences  due  to  its  neighbors  is  summed  to  obtain  the  resultant  force  for  that 
particle. 

The  range  of  the  attractive  force  for  Cu— Cu  is  assumed  to  be  2.4  lattice 
units.  This  range  corresponds  to  a  point  between  the  second  and  third  nearest 
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neighbors.  This  range  cut-off  has  been  investigated  and  found  to  be  reasonable  for 
fee  targets  [ref.  56].  The  differences  in  potential  beyond  fhe  second  nearest 
neighbors  were  found  to  be  minor. 


2.  Genera]  Forms 

The  atom— atom  potentials  are  composite  functions  consisting  of  a 
repulsive  section,  the  wall,  joined  smoothly  to  an  attractive  section,  the  well.  The 
joining  of  the  two  functions  is  accomplished  by  a  cubic  spline  function.  The 
atom— ion  potential  is  a  strictly  repulsive  wall,  which  is  much  steeper  than  that  of 
the  atom— atom. 

The  atom— atom  repulsive  wall  is  a  Born— Mayer  potential  function  given 
by  [ref.  57]. 

V(r)  =  A  e~Br 

where  A  is  the  atom's  hardness  and  R  is  the  atom's  size. 

The  atom-ion  repulsive  wall  is  a  Moliere  potential  function  given  by, 

V(r)  =  [(Z,Z2  e2/a)/(r/a)]ip(r/a) 

where 


VK r/a)  =  [0.35e~°'3r/a  +  0  55e  L2r/a  +  0.1e“6'0r/a] 


a  =  0.8853a, 


•  a .  =  Bohr  radius 

O 


.r~zi+r^2. 


-2/3 


•  Zi  =  atomic  number  of  atom 

•  Z2  =  atomic  number  of  ion. 
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The  general  form  used  for  the  well  is  a  Morse  potential  of  the  form, 
V(r)  =  DJ  e-<2“(r-re)>  -2e-<Q(r-r«>] 

where 


•  De  =  well  depth 

•  re  =  equilibrium  separation 


•  a  —  well  width  control. 

1.  The  Functions  Used  in  the  Simulation 

The  atom— atom  potential  is  a  repulsive  Born— Mayer  joined  to  an 
attractive  Morse  potential  by  a  cubic  spline. 

The  composite  pair  potential  function  is  given  bv, 


Vjj.  =  Co+Cir  +  C2r  +  C3r3 

Vjj  =  DJ  e-<2a(i— re)>  _  2e-<o(r-re)>] 


r  <  Ra 

Ra  <  r  <  Rb 
Rb  <  r  <  Rc 
Rc  <  r 


where  is  the  composite  potential  for  the  ilh  and  jth  atom  separated  by  a 
distance  r. 

The  cubic  spline  joins  the  Born-Mayer  and  Morse  at  points  Ra  and  Rb. 
The  attractive  Morse  function  is  truncated  at  Rc,  between  the  second  and  third 
nearest  neighbors. 

The  atom— ion  potential  is  obtained  with  the  purely  repulsive  Moliere 


function, 


vi  =  [(Z1Z2  e2/a)  /  (r/a)]  <^( r/a)  r  <  Ra 

Vj  =  0  r  >  Ra 

where  V;  is  the  potential  of  the  ith  atom  and  the  ion. 
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The  potential  function  parameters  used  in  the  simulation  are  found  in 
Table  2  of  Appendix  B.  The  cubic  spline  parameters  are  found  in  Table  3  of 
Appendix  B.  The  actual  Ar-Cu  and  Cu-Cu  potential  functions  are  plotted  in 
Figures  7  and  8  of  Appendix  A. 

E.  ERROR  ANALYSIS 

1 .  The  Experimental  Error 

Since  the  system  is  non  dissipative,  energy  is  conserved  and  the  energy- 
calculation  results  may  be  used  as  a  possible  check  of  numerical  error.  In  sputtering 
simulations,  a  small  percent  of  energy  error  (up  to  3%)  is  considered  acceptable. 
The  justification  is  that  studies  have  shown  that  a  small  change  in  the  initial  ion 
energy  does  not  significantly  affect  the  sputtering  yield.  In  the  simulation  of  a 
liquid,  however,  the  conservation  of  energy  is  an  important  factor  which  should  be 
maintained. 

The  principal  reason  for  our  concern  with  energy  conservation  is  that  an 
error  in  energy  is  related  to  an  error  in  temperature  by, 

energy  error  =  3/2Nk  x  temperature  error 
where  N  is  the  number  of  particles  and  k  is  the  Boltzmann  constant. 

This  equation  says  that  a  1  eV  uncertainty  in  energy  corresponds  to  a 
•  maximum  uncertainty  in  temperature  of  5.35  K.  This  means  that  a  1%  error  in  the 
total  energy  of  a  liquid  target  ,  say  2461  eV,  equilibrated  to  a  temperature  of  1575  K 
would  correspond  to  an  uncertainty  in  temperature  of  ±  65.95  K.  Such  an 
uncertainty  in  temperature  would  be  critical  if  the  temperature  of  the  liquid  is  close 
to  either  the  melting  or  boiling  points.  The  largest  energy  loss  in  the  simulation 
was  6.45  eV  with  a  corresponding  temperature  uncertainty  of  34.51  K. 
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2.  Establishing  the  Fit  of  Curves 

According  to  Maxwell's  equipartition  of  energy,  the  energy  added  to  a 
number  of  molecules  will  be  equally  partitioned  between  the  potential  energy  and 
the  components  of  kinetic  energy  as  the  ensemble  comes  to  thermal  equilibrium. 
The  variation  of  total  kinetic  energy  with  time  is,  therefore,  a  good  indicator  of 
thermal  equilibrium.  In  our  simulation,  when  the  total  kinetic  energy  approaches  a 
value  equal  to  half  of  the  total  added  energy  and  when  the  components  of  kinetic 
energy  become  more  or  less  equal,  the  group  of  atoms  will  be  considered  to  be  in 
thermal  equilibrium.  How  close  these  relationships  come  to  their  theoretically 
expected  values  will  be  the  principal  measure  of  target  equilibration. 

The  third  type  of  relationship  generated  from  the  results  of  the  liquid 
simulation  programs  is  the  radial  distribution  of  the  atoms.  This  has  been 
determined  experimentally  [ref.  46].  The  calculated  radial  distribution  function  will 
be  compared  to  the  experimental  radial  distribution  function  by  the  general  shape 
and  by  the  coincidence  of  their  peaks  and  valleys. 

The  velocity  distribution  can  be  compared  to  the  theoretical  expected 
value  using  statistical  analysis.  The  expected  distribution  of  the  atom  velocities  is 
the  Maxwell-Boltzmann  distribution.  \2  tests  will  be  used  to  determine  how  well 
the  calculated  velocity  distributions  fit  the  expected  theoretical  values  of  the 
Maxwell-Boltzmann  distribution  function.  The  following  equation  is  used  to  obtain 
the  x2  values, 

e=i/jY(°>  ~  Et>2 
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where 


•  d  =  numberof  degrees  of  freedom 

•  k  =  1 • • -n,  individual  number  of  bins 

•  0^  =  observed  values 

•  Ek  =  expected  values. 

Generally,  if  x2 » U  the  observed  results  do  not  fit  the  assumed 
distribution  very  well.  The  fit  is  generally  considered  acceptable  for  a  \2  <  1. 

A  more  qualitative  measurement  of  the  fit  is  made  by  finding  the 
percentage  probability  Pd(x2  >  X'2)  of  obtaining  a  value  of  \ 2  >  Aq-  Where  x2  is  the 
actual  value  of  x2  obtained  from  experimental  measurements.  These  probabilities 
are  found  in  Table  D,  page  251  of  reference  58.  If  the  probability  Pd(\/2  >  X2)  <  5% 
then  the  fit  is  considered  to  be  acceptable.  This  is  the  criteria  applied  to  our 
results. 

In  our  simulation,  the  number  of  bins  used  for  the  velocity  distribution  is 
30  and  the  degrees  of  freedom  is  d  =  29.  Using  Table  D  of  reference  58,  the 
maximum  x-2  corresponding  to  a  Pd(x2^  Aq)  =  5 %  for  29  degrees  of  freedom  is 
1.493.  Therefore,  if  the  x2  of  any  of  our  velocity  distributions  is  under  1.493,  then 
the  fit  will  be  considered  to  be  consistent  with  the  expected  values. 
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IV.  RESULTS  AND  DISCUSSION 


Eight  different  targets  were  generated  using  various  boundary  conditions, 
warming  methods  and  nearest  neighbor  list  update  schemes.  The  target  with  the 
best  liquid-like  characteristics  was  chosen  and  subjected  to  a  295  trajectory 
sputtering  run. 

In  order  to  facilitate  reference  to  the  various  simulation  schemes  and  their 
results  the  following  names  are  defined, 


•  QLP,  Liquid  simulation  using  the  atom  displacement  warming  scheme  and 
unrestrained  boundaries. 

•  QLPBC,  Liquid  simulation  using  the  atom  displacement  warming  scheme 
and  reflective  boundary  conditions. 

•  QLP-REV1,  Liquid  simulation  using  the  atom  displacement  warming 
scheme,  unrestrained  boundaries  and  a  new  neighbor-list  update  every 
time  step. 

•  QLPBC-REV1,  Liquid  simulation  using  the  atom  displacement  warming 
scheme,  reflective  boundary  conditions  and  a.  new  neighbor— list  update 
every  time  step. 

•  QLV,  Liquid  simulation  using  the  impulse  warming  scheme  and  unrestrained 
boundaries. 

•  QLVBC,  Liquid  simulation  using  the  impulse  warming  scheme 
and  reflective  boundary  conditions. 

•  QLV-REV1,  Liquid  simulation  using  the  impulse  warming  scheme, 
unrestrained  boundaries  and  a  new  neighbor— list  update  every  time  step. 

•  QLVBC— REV1,  Liquid  simulation  using  the  impulse  warming  scheme, 
reflective  boundary  conditions  and  a  new  neighbor-list  update 

every  time  step. 
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BC  means  that  reflective  boundary  conditions  are  applied  and  REV1  stands  for 
the  revision  which  includes  the  update  of  the  nearest  neighbor  list  every  time  step. 

A.  ENERGY  CONSERVATION 

The  calculated  total  energy  losses  for  all  eight  runs  range  from  0.001%  to  0.25%. 
The  equivalent  energy  uncertainty  for  this  range  is  0.03  eV  to  6.45  eV.  Although  all 
eight  results  have  acceptable  energy  losses,  there  is  a  clear  division  in  the  energy 
conservation  efficiency  betweeen  the  runs  that  update  the  nearest  neighbor  list  every 
time  step  and  those  that  do  not.  The  simulation  runs  which  update  the  nearest 
neighbor  list  each  time  step  give  the  poorest  energy  conservation.  This  is  to  be 
expected  because  as  the  simulation  proceeds,  atoms  find  new  neighbors  with  which 
they  can  interact  and  this  requires  additional  calculations  that  add  to  the  total 
error. 

The  best  energy  conservation  was  obtained  with  QLPBC.  A  summary  of  the 
experimental  error  for  all  eight  targets  is  found  in  Table  4  of  Appendix  B. 

B.  DENSITY  RESULTS 

The  density  was  calculated  for  each  of  the  unrestrained  targets  (without 
reflective  boundary  conditions)  after  equilibration  was  reached.  These  were  found 
to  be  slightly  lower  than  those  corresponding  to  liquid  copper  at  the  equilibration 
temperature. 

The  most  accurate  density  was  obtained  with  QLV-REV1.  This  simulation 
yielded  a  liquid  density  of  7.664  gm/cm3,  within  1.8%  of  the  expected  value  for  a 
liquid  temperature  of  1690  K.  Table  5  of  Appendix  B  shows  the  resulting  densities 
for  all  runs  with  unrestricted  boundaries  and  their  expected  theoretical  values. 
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C.  TARGET  EQUILIBRATION  AND  FINAL  TEMPERATURES 

The  plots  of  total  kinetic  energy  versus  elapsed  time  (Figs.  9—16,  Appendix  A) 
and  those  of  the  components  of  kinetic  energy  versus  elapsed  time 
(Figs.  17—24,  Appendix  A)  show  that  all  targets  reach  equilibrium  after 
approximately  400  femptoseconds.  No  significant  change  in  the  kinetic  energy  or  in 
the  components  of  kinetic  energy  is  observed  after  this  time.  This  equilibration 
time  seems  to  be  consistent  with  the  simulation  results  of  Lo  et  al  [ref.  42],  The 
simulations  were  computed  for  at  least  twice  the  equilibration  time. 

These  kinetic  energy  plots  also  show  that  equipartition  of  energy  holds  very  well 
for  all  eight  target  generation  schemes.  In  particular,  QLPBC  produced  a  final 
kinetic  energy  of  273  ±  17  eV,  less  than  1  eV  difference  from  the  theoretical  expected 
value  of  273.071.  The  error  for  the  total  kinetic  energies  is  taken  from  the 
maximum  amplitude  of  the  kinetic  energy  oscillations  about  an  average  value 
following  the  equilibration  point.  Table  6  of  Appendix  B  summarizes  the  final 
kinetic  energy  and  temperature  results  for  all  eight  targets. 

QLV— REVl  and  QLVBC— REV1  produce  the  smallest  kinetic  energy 
oscillations  at  equilibrium,  ±4  eV.  This  corresponds  to  a  temperature  uncertainty  of 
±32  K.  The  final  temperatures  for  these  two  targets  are  1690  ±  32  I\  and  1700  ±  32 
respectively.  The  reason  for  such  oscillations  in  kinetic  energy  after  equilibration  is 
reached  is  that  the  size  of  the  ensemble  is  small  compared  to  a  real  liquid  system. 
Real  systems  have  many  billions  of  particles  and  the  total  kinetic  energy 
fluctuations  can  average  out. 
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D.  RADIAL  DISTRIBUTIONS 

The  radial  distribution  for  an  unwarmed  solid  copper  crystal  is  given  by 
Figure  25  of  Appendix  A.  A  comparison  of  this  figure  with  the  various  radial 
distributions  produced  by  the  liquid  simulations  (Figs.  26—33  of  Appendix  A) 
clearly  contrasts  the  structure  differences  between  a  solid  and  a  liquid.  These  liquid 
radial  distributions  have  the  distinctive  wide  and  smooth,  peaks  and  valleys 
characteristic  of  a  liquid.  The  crystal  structure  is  completely  lost  in  all  samples. 

The  peaks  and  valleys  of  all  the  resulting  radial  distributions  are  compared 
against  each  other  for  consistancy,  and  against  neutron  diffraction  experimental 
results  [ref.  46]  for  validation,  in  Table  7  of  Appendix  B.  All  the  numbers  for  this 
table  are  taken  directly  from  the  radial  distributions  plots,  Figures  26-33  of 
Appendix  A.  It  is  important  to  note  that  many  of  the  expected  valleys  could  not  be 
positively  identified  and  were  therefore  left  blank. 

The  best  (most  liquid-like)  radial  distributions  are  those  obtained  with 
QLP-REV1.  QLPBC— REV1,  QLV-REVl,  and  QLVBC-REVl.  These  radial 
distributions  are  in  excellent  agreement  with  the  neutron  diffraction  experimental 
results  [ref.  46].  The  very  best  radial  distribution  of  all  is  that  of  QLVBC-REVl. 

E.  VELOCITY  DISTRIBUTIONS 

The  resulting  velocity  distribution  and  the  expected  Maxwellian  distribution  for 
each  target  are  plotted  together  in  Figures  34—41  of  Appendix  A.  These  figures 
show  that  the  velocity  distribution  of  all  eight  liquid  targets  generally  match  the 
Maxwellian  distribution.  The  reduced  x2  tests  and  their  probabilities  Pd(x2^X2)  are 
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found  in  Table  8  of  Appendix  B.  This  table  shows  that  the  most  Maxwellian— like 
velocity  distribution  is  obtained  with  QLP.  Only  three  targets  QLP,  QLV  and 
QLVBC— REV1  have  velocity  distributions  with  a  Pd(\2>\2)  >  5%.  This  is  the 
standard  rejection  criteria.  The  smallest  Pd(x2^\'2)  was  0.5%  with  most  falling 
around  1.5%.  In  general,  the  velocity  distributions  appear  to  be  consistent  with  the 
Maxwellian  distribution. 

F.  SPUTTERING  RESULTS 

The  QLVBC— REV;1  target  was  chosen  for  the  sputtering  run  because  it  has  the 
best  overall  liquid— like  characteristics.  295  trajectories  were  run  with  1.0  KeV 
argon  ions.  The  liquid  target  atom  velocities  were  zeroed  to  save  both 
computational  costs  and  real  time  running  requirements.  A  sputtering  run  with  live 
atoms  (all  atoms  initially  moving)  takes  on  average  about  five  times  longer  than  a 
static  target. 

The  sputtering  program  had  on  average  a  3%  energy  loss  per  trajectory. 
Although  this  is  considered  acceptable  for  a  sputtering  simulation,  it  was  later  found 
that  the  program  overestimated  the  amount  of  energy  that  left  the  target  with  every 
ejected  atom,  but  did  not  affect  the  forces.  Thus  3%  is  very  much  an  upper  bound. 

The  final  average  yield  for  the  295  trajectories  is  8.28  atoms  per  incident  ion. 
This  is  40%  higher  than  the  typical  yield  expected  from  a  solid  target  on  the  010 
face  and  at  room  temperature.  This  increase  in  yield  is  in  qualitative  agreement 
with  experimental  results  [refs.  33-35]  and  with  another  computer  simulation 
[ref.  43].  Table  9  in  Appendix  B  shows  the  general  consistancy  of  these  results. 
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V.  CONCLUSIONS 


Results  from  the  various  liquid  target  generation  schemes  have  shown  that  the 
warming  method  used  to  energize  the  individual  target  atoms  does  not  affect  the 
final  liquid— like  characteristics  of  the  resulting  target.  That  is,  results  of  liquid 
targets  warmed  by  impulse  were  essentially  indistinguishable  in  radial  distribution 
and  in  velocity  distribution  from  those  warmed  by  displacement.  Both  warming 
methods  are  equally  capable  of  generating  liquid  targets. 

The  best  results  were  obtained  with  the  application  of  reflection  boundary 
conditions  in  conjunction  with  a  scheme  which  updated  the  nearest  neighbor  list 
every  time  step.  QLPBC— REVl  and  QLVBC— REVl  produced  good  liquid  targets. 
The  liquid  target  generated  by  QLVBC-REV1  is  found  to  be  the  best  of  the  two. 
Figure  43  shows  the  actual  liquid  structure  generated  by  QLVBC-REVl.  The 
straight  lines  in  this  figure  represent  the  outline  of  the  crystal  before  warming  and 
equilibration. 

Results  indicate  that  the  motion  of  a  typical  liquid  atom  in  a  time  span 
equivalent  to  at  least  twice  the  equilibration  time  (approximately  800  fs)  is  such 
that  it  seldom  leaves  its  original  atom  neighborhood.  Figure  42  in  Appendix  B 
illustrates  this  behavior. 

The  sputtering  of  the  best  liquid  target  with  a  1  KeV  incident  argon  ion 
resulted  in  yields  which  were  40%  higher  than  expected  by  solid  targets.  This  is  in 
qualitative  agreement  with  the  results  of  other  experimental  efforts  [refs.  33-35] 
found  in  Table  9  of  Appendix  B.  Figures  44  and  45  in  Appendix  A  show  the  yield 
per  incident  particle  distribution  for  both  a  crystal  Cu(010)  target  and  a  Liquid  Cu 
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target  respectively.  The  large  number  of  trajectories  that  have  no  yield  at  all  for 
the  crystal  target  may  be  attributed  to  ion  channeling. 

The  spot  pattern  (spatial  distribution  of  ejected  material)  of  the  liquid  target, 
Figure  46,  shows  no  signs  of  an  ordered  structure.  The  spot  pattern  of  a  Cu(010) 
crystal,  Figure  47,  is  shown  to  contrast  the  difference  in  spot  pattern  between  a  well 
ordered  structure  and  an  amorphous  ensemble. 

In  summary,  the  liquid  target  generation  methods  QLVBC-REVl  and 
QLPBC— REVl  both  produced  reasonable  liquid  targets.  The  frequent  update  of  the 
nearest  neighbor  list  as  well  as  the  application  of  boundary  conditions  was  found  to 
he  critical  in  making  a  reasonable  liquid  structure.  The  sputtering  results  are  in 
general  agreement  with  published  data. 
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VI.  RECOMMENDATIONS 


The  primary  purpose  of  this  investigation  was  to  devise  a  method  that  would 
produce  a  reasonable  liquid  target.  This  was  achieved. 

The  secondary  objective  wras  to  conduct  a  preliminary  study  of  liquid  sputtering 
using  the  best  liquid  target.  The  average  sputtering  yield  was  obtained  from  a  set 
of  295  trajectories  using  only  one  liquid  target.  Other  targets  should  also  be  used  in 
sputtering  runs  in  order  to  obtain  additional  evidence  for  the  validity  of  the  liquid 
targets.  It  is  possible  that  by  using  the  same  liquid  target  for  every  iiajectory  we 
may  be  introducing  a  systematic  bias  in  the  sputtering  yield.  It  is  reasonable  to 
assume,  for  example,  that  our  yields  are  closer  to  that  of  an  amorphous  solid  than 
that  of  a  liquid  because  the  target  has  constant  structure. 

A  scheme  should  be  devised  by  which  several  equivalent  liquid  targets  an  used 
in  turn  (cycled)  with  each  trajectory  of  a  full  sputtering  run.  Another  alternative 
would  be  to  complete  several  sets  of  trajectories  on  sufficiently  separated  areas  of 
the  surface. 

At  least  one  liquid  target  should  be  subjected  to  various  ion  energies  to 
determine  if  there  is  an  energy  dependent  yield  profile  as  suggested  by  the  results  of 
reference  33.  The  effect  on  yield  by  varying  the  angle  of  incidence  should  also  be 
investigated. 


43 


APPENDIX  A  -  FIGURES 


Figure  1.  The  Rebound  Sputtering  Process 


Figure  2.  The  Recoil  Sputtering  Process 
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Figure  3.  The  Reflection  Sputtering  Process 


Figure  4.  The  Direct  or  Deflected  Recoil  Process 


Figure  5.  The  Collision  Cascade 
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Figure  6.  The  Target-Ion  Alignment 
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Figure  7.  Ar— Cu  Potential  Function. 
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Figure  8.  Cu-Cu  Potential  Function. 
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Figure  10.  Kinetic  Energy  Versus  Elapsed  Time,  QLPBC  Output. 
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Figure  11.  Kinetic  Energy  Versus  Elapsed  Time,  QLP-REVl  Output. 
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Figure  12.  Kinetic  Energy  Versus  Elapsed  Time,  QLPBC— REV1  Output. 
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Figure  13.  Kinetic  Energy  Versus  Elapsed  Time,  QLV  Output,. 
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Figure  15.  Kinetic  Energy  Versus  Elapsed  Time,  QLV— REVl  Output. 
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Figure  16.  Kinetic  Energy  Versus  Elapsed  Time,  QLVBC— REV1  Output. 
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Figure  18. 


Components  of  KE  Versus 


Elapsed  Time,  QLPBC  Output. 
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Figure  23.  Components  of  KE  Versus  Elapsed  Time,  QLV— REVl  Output. 
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Figure  24.  Components  of  KE  Versus  Elapsed  Time,  QLVBC— REV1  Output. 
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Figure  25.  Radial  Distribution,  Unwarmed  Crystal. 


65 


RADIAL  SEPARATION  (  angstroms  ) 


RADIAL  DISTRIBUTION 


I  I  I  I  I  I 


n 

to 

l£> 

CN 

CN 

CO 

ID 

CN 

<— 

CO 

•*- 

CN 

o 

d 

rsi 

CN 

Csj 

csj 

d 

T - 

r— 

d 

o 

q 

q 

q 

d 

O 

d 

o 

o 

d 

d 

d 

d 

o’ 

d 

d 

Ainiave  oud 


Figure  26.  Radial  Distribution,  QLP  Output. 
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Figure  27.  Radial  Distribution,  QLPBC  Output. 
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Figure  30.  Radial  Distribution,  QLV  Output. 
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Figure  31.  Radial  Distribution,  QLVBC  Output. 
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Figure  32.  Radial  Distribution,  QLV-REVI  Output. 
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Figure  33.  Radial  Distribution,  QLVBC-REV1  Output. 
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Figure  34.  Velocity  Distribution.  QLP  Output. 
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Figure  36.  Velocity  Distribution,  QLP— REV1  Output. 
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Figure  37.  Velocity  Distribution,  QLPBC— REV1  Output. 
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Figure  38.  Velocity  Distribution,  QLV  Output. 
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Figure  39.  Velocity  Distribution,  QLVBC  Output. 
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Figure  45.  Ejected  Atoms  per  Single  Ion  for  Liquid  Target  QLVBC-REV1 
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Figure  47.  Spot  Pattern  for  Solid  Target,  Cu(010)  Crystal 
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APPENDIX  B  -  TABLES 


TABLE  1 

PHYSICAL  DATA  FOR  COPPER 

Atomic  Number  (Z) 

29 

Atomic  Weight 

63.540  amu 

Crystal  Type 

FCC 

Lattice  Constant 

3.6150  angstroms 

Re 

1.8075  angstroms 

Melting  Point 

1083.4  ±  0  2  oc 

Boiling  Point 

2567  oc 

Density  (solid) 

8.96  gm/cm3 

Note:  Data  for  Table  1  is  derived  from  reference  59. 
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TABLE  2 

POTENTIAL  FUNCTION 

PARAMETERS 

Parameter 

Cu-Cu 

Cu-Ar 

A  (KeV) 

22.564 

B  (A'1) 

-5.088 

De  (eV) 

0.431 

Re  (A) 

2.628 

a  (A-*) 

1.405 

Ra  (LU) 

0.83 

1.41 

Rb  (LU) 

1.10 

1.41 

Rc  (LU) 

2.40 

1.41 

Zi  (amti) 

29 

18 

Zo  (amu) 

29 

29 

K 

0.0 

0.0920 

TABLE  3 

CUBIC  SPLINE  PARAMETERS 

j  Parameter 

Cu-Cu 

Co 

5.876xl02 

c, 

-1.594x103 

c2 

1.450xl03 

c3 

^.423xl02 
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TABLE  4  [ 

TOTAL  ENERGY 

LOSSES 

Target 

Percent  Error(%) 

Energy  Uncertainty^  V)  | 

QLP 

0.001 

0.04 

QLPBC 

<  0.001 

0.03 

QLP-REV1 

0.251 

6.44 

QLPBC-REVl 

0.238 

6.45 

QLV 

0.01 

0.30 

QLVBC 

0.01 

0.29 

QLV-REV1 

0.23 

6.22 

QLVBC-REV1 

0.24 

4.20 

TABLE  5 

DENSITY  RESULTS  (TARGETS  WITH  I 

JNRESTRICTED  BOUNDARIES)  ( 

Density  (gm/cm3) 

Target 

Expected 

Result 

QLP 

7.906 

6.646 

QLP-REV1 

7.906 

6.617 

QLV 

7.806 

6.057 

QLV-REYT 

— 

7.806 

7.364 
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TABLE  6 


KINETIC  ENERGY  AND  TEMPERATURE  RESULTS 


Target 

KE  (eV) 

Temperature  (I<)  j 

Expected 

— 

Expected 

Result 

QLP 

273.071 

265  ±  6 

1462.043 

1420  ±  32 

QLPBC 

273.071 

273  ±  6 

1462.043 

1460  ±  32 

QLP-REV1 

273.071 

280  ±  17 

1462.043 

1520  ±  91 

QLPBC-REV1 

273.071 

290  ±  8 

1462.043 

1560  ±  43 

QLV 

296.458 

288  ±  9 

1587.258 

1540  ±  48 

QLVBC 

296.458 

299  ±  12 

1587.258 

1600  ±  64 

QLV-REVl 

296.458 

315  ±  4 

1587.258 

1690  ±  21 

QLVBC-REVl 

296.458 

280  ±  4 

1587.258 

1700  ±  21 
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TABLE  7 


RADIAL  DISTRIBUTIONS,  PEAKS  AND  VALLEYS 


Target 


1st 

Peak 

(LU) 


1st 

2nd 

pppM 

3rd 

3rd 

Valley 

Peak 

.  I 

Peak 

Valley 

(LU) 

(LU) 

■0 BM 

(LU) 

(LU) 

QLP 

QLPBC 

QLP-REVl 

QLPBC-REVI 

QLV 

QLVBC 

QLV-REVl 

QLVBC-REVl 


2.5±.l 

2.5±.l 

2.4±.l 

2.4±.l 

2.5±.l 

2.5±.l 

2.4±.l 

2.4±.l 


3.3±.l 


3.3±.l 

3.3±.l 


4.3±.l 

4.7±.l 

4.3±.i 

4.4±.l 

4.3±.l 

4.3±.l 

4.3±.l 

4.3±.l 


5.6±.l 


6.7±.l 
6.7±.l 
6.6±.l 
6.6±.l 
6.7±.  1 
6.7±.l 
6.6±.l 
6.5±.l 


7.7±.l 

7.7±.l 


NEUTRON 
DIFF.  DATA 
[ref.  46] 


2.5±.I 


3.5±.l 


4.7±.l 


5.6±.l 


6.8±.l 


7.9±.l 


TABLE  S  | 

VELOCITY  DIST. 

V2  TESTS  TO  MAXWELLIAN  DIST. 

Target 

\2 

Pd(\2>.\2o) 

QLP 

0.9 

61.5 

QLPBC 

1.7 

1.6 

QLP-REY 1 

1.7 

1.6 

QLPBC-REY1 

1.9 

0.5 

j  QLV 

1.5 

5.4 

QLVBC 

1.7 

1.6 

QLY-REYl 

1.5 

5.4 

QLYBC-REY1 

1.5 

1.5 

TABLE  9 

PERCENT  INCREASE  IN  YIELD  1 

OF  LIQUID  SPUTTERING  OYER  SOLID  SPUTTERING  j 

Source 

Type  of 
Study 

Target 

Ion 

Energy 

(KeV) 

— - - 1 

Percent  Yield 

Increase  (%) 

Ref.  33 

experiment 

liquid  Sn 

Argon 

40 

-6 

Ref.  34 

experiment 

liquid  Sn 

Argon 

0.2 

0.375 

50 

15 

Ref.  35 

experiment 

liquid  In 

Argon 

0.107 

10 

Ref.  43 

simulation 

liquid  Cu 

Argon 

5.0 

60 

This 

Study 

simulation 

liquid  Cu 

Argon 

1.0 

40 
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